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ABSTRACT 

The gravitational instability of a dust layer is one of the scenarios for planetesi- 
mal formation. If the density of a dust layer becomes sufficiently high as a result of the 
sedimentation of dust grains toward the midplane of a protoplanetary disk, the layer be- 
comes gravitationally unstable and spontaneously fragments into planetesimals. Using 
a shearing box method, we performed local A^-body simulations of gravitational insta- 
bility of a dust layer and subsequent coagulation without gas and investigated the basic 
formation process of planetesimals. In this paper, we adopted the accretion model as a 
collision model. A gravitationally bound pair of particles is replaced by a single particle 
with the total mass of the pair. This accretion model enables us to perform long-term 
and large-scale calculations. We confirmed that the formation process of planetesimals 
is the same as that in the previous paper with the rubble pile models. The formation 
process is divided into three stages: the formation of non-axisymmetric structures, the 
creation of planetesimal seeds, and their collisional growth. We investigated the depen- 
dence of the planetesimal mass on the simulation domain size. We found that the mean 

3 /2 

mass of planetesimals formed in simulations is proportional to Ly , where Ly is the size 
of the computational domain in the direction of rotation. However, the mean mass of 
planetesimals is independent of Lx, where is the size of the computational domain 
in the radial direction if is sufficiently large. We presented the estimation formula 
of the planetesimal mass taking into account the simulation domain size. 

Subject headings: gravitation, instabilities, methods:n-body simulation, planets and 
satellites:formation 
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Introduction 



According to the standard theory of planet formation, planets form from planetesimals, which 
are kilometer-sized solid bodies. However, the process of planetesimal formation is still controver- 
sial. The i nward migration of meter-sized bodies du e to gas drag is very rapid; its timescale is 
only 10^ yr ( Adachi et al.lll976l : IWeidenschillinglll977l ). The growth of dust to planetesimals due to 
particle adhering seems difficult in the standard disk model. 



The gravitational in stability scenario is an alternative scenario for this stage (|Safronovlll969l : 
Goldreich &: Wardll 19731 ). Dust particles settle into the midplane of a protoplanetary disk owing to 
the gravitational force from the central star if the gas flow is laminar. As the sedimentation of dust 
particles proceeds, the density at the midplane exceeds the critical density and the dust layer finally 
becomes gravitationally unstable. The gravi tationally unstable dust l ayer rapidly cohapses by its 



self-g ravity, forming km-sized planetesimals (j Goldreich &: Ward 



1973 



Coradini et al 



1981 



Sekiva 



19831 ). The rapid migration of meter-sized bodies can be avoided because the timescale of the grav- 
itational instability is only on the order of a Kepler period. However, the gravitational instability 
scenario has a critical issue. As the sedimentation of dust grains toward the midplane proceeds, 
the vertical velocity shear increases and gives rise to the Kelvin-Helmholtz instability, which makes 



ies have been carried out on this issue ( Weidenschilline 


1980: 


Cuzzi et al. 


1993; 


Chamonev et al. 


1995 




Weidenschilhnel Il995l: 


Sekiva 


1998 




Dobrovolskis et al. 


1999: 


Sekiva &: Ishitsu 


2000. 


2001; 


Ishitsu &: Sekiva 


2002. 


2003: 


Michikoshi & Inutsuka 


2006 


). However, problems concerning the grav- 



itational instability still remain unsettled. 

Although many studies have been conducted on the condition of gravitational instability or its 
linear regime, there has been little study of the nonlinear regime of gravitational instability. Here 
we concentr a te oii the formation process of planetesimals assuming the gravitational instability. 
Tanga et al.l (|200J) performed A^-body simulations of a gravitationally unstable particle disk. In 
their simulation, the optical depth is very small, and the particle disk is unstable, even initially. If we 
consider the formation of planetesimals through gravitational instability at ao — lAU, the optical 
depth is high, and the initial velocity dispersion must be large o wing to turbulent flow d riven by 

In our 



1991 



Johansen et al 



the Kelvin-Helmholtz instability or magnetorotational instability (iBalbus Sz Hawley 
analys is, we will take a close look at initially stable disks with large optical depth. 
([20071) performed simulations using a code that solves the magnetohydrodynamic equations with a 
three-dimensional grid and includes particles with self-gravity. They mapped the particle density 
on the grid and solved the Poisson equation using a fast Fourier transform method. They included 
collisions among particles as damping the velocity dispersion of the particles within each grid. They 
found that particles collapse in an overdense region in the midplane and the gravitationally bound 
objects form with masses comparable to dwarf planets. 



Michikoshi et al.l ()2007l ) performed a set of simulations in the self-gravitating collision-dominated 



particle disks without gas components (hereafter Paper I). The point-to-point Newtonian mutual 
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interaction was calculated directly. They adopted the hard and soft sphere models as collision 
models, and examined parameter dependencies on the size of computational domain, restitution 
coefficient, optical depth. Hill radius of particles, and the d uration of collision. In the hard sphere 



model, penetration between particles is not possible (e.g. , IWisdom &: Tremaind 119881 : 1 Said Il991 



Richar dsoni 1 1 994 : iDaisaka Sz Idalll999l : iDaisaka et al.ll200ll ). The yelocity changes in a n instant a t 
collisions. In the soft sphere model, particles can overlap if they are in contact (e.g., ISaldll995l ). 
They found that the formation process of planetesimals is qualitatively independent of simulation 
parameters if initial Toomre's Q value, Qinit > 2. The forma tion proce ss is divided into three 
stages: the formation of non-axisymmetric wake- like structures (|Saldll995l ). the creation of aggre- 
gates, and the rapid collisional growth of the aggregates. The mass of the largest aggregates is 
larger than the ma.ss pre dicted by t he linear perturbation theory (hereafter linear mass Miinear) 
(jGoldreich WardI Il973l ) . However iMichikoshi et al.l (|2007l ) could not find the saturation of the 
growth of the planetesimals in Paper I. The mass of the planetesimal formed in the simulation 
depends on the size of the computational domain. Almost all mass in the computational domain 
is absorbed by a single planetesimal, and thus the mass of the planetesimal is determined by the 
total mass in the computational domain. 

In this paper, we use the alternative model of collisions, 'accretion model'. In physical terms, 
the accretion model corresponds to the more dissipative model than the hard and soft sphere 
models. In the accretion model, when two particles collide, if the binding condition is satisfied, the 
two particles are merged into one particle. Efficient energy dissipation occurs when two particles 
merge into one particle. In the soft and hard sphere models, since large aggregates consist of a lot of 
small particles and the number of particles does not decrease, the calculation is time-consuming. If 
we are not interested in internal states of planetesimals such as rotation or internal density profile, 
we can treat a large aggregate as one large particle. We expect that the final mass or spatial 
distribution of planetesimals can be estimated by the accretion model. The accretion model has an 
advantage from the viewpoint of calculation. The number of particles decreases as the calculation 
proceeds; thereby this model enables us to perform simulations of larger numbers of particles than 
that in the previous work. Taking this advantage, we can examine the large computational domains. 

In this study, we neglect the effect of gas for the sake of simplicity. As many authors in- 
yestigated, the effect of gas is very i mpor tant. The solid particles drift radially due to gas dra; 



(Adachi et al 



19761 : IWeidenschilling] 119771 ) , gas friction helps gravitatio nal collapse (IWard 



I993I : barge k Sommerial Il995l : 



Fromang &: Nelson 



solid particles form d ue to streaming instability (jYoudin &: Goodman 



2006 



Johansen et al 



YoudiJboOSallbl'l. gas turbulence prevents and helps concen t ration of solids (^Weidenschilling &: C 



1971 



2005 



2006al). and clumps of 



Johansen et al.ll2006d : 



Johansen et al.l 120071 ). But once gravitational instability happens and large aggregates form, as a 
first step, we may be able to neglect gas drag because they are so large that they decouple from 
gas. Therefore our gas-free model may be applicable to the stages of gravitational instability and 
subsequent collisional growth. 



In Section [21 we explain the simulation method. Our results are presented in Section [3j We 
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compare the accretion model with the soft and hard sphere models and investigate the models of 
large computational domain. In Section 31 we summarize the results. 



2. Methods of Calculation 



The method of calculation is the same as that used in Paper I ex cept for the collision model 



which is based on the method of simulation of planet ary rings (e.g., IWisdom &: Tremaing 
Salolll99ll : lOaisaka Idalll999l : lohtsuki k Emorilboool ). 



198S 



We introduce rotating Cartesian coordinates (x, y, z) called Hill coordinates, the x-axis pointing 
radially outward, the y-axis pointing in the direction of rotation, and the z-axis normal to the 
equatorial plane. The origin of the coordinates moves on a circular orbit with the semi-major axis 
ao and the Keplerian angular velocity i7o- Assuming [zjl <C ao, where {xj,yj,Zj) is the 

position of jth particle, we can write the equation of motion in non-dimensional forms independent 
of the mass and the semi-major axis oq, if we scale the time by , the length by Hill radius hao, 



and the mass by h?M^, where h = (2m p/3M. 
is the mass of the central star (e.g. 



a/3 



Petit fc Henon 



is the characteristic mass of particles, and 



1986 



Nakazawa fc Idalll988l l: 
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where a tilde denotes corresponding non-dimensional variables, rij is the distance between particles 
i and j, and is the number of particles. 



We use the periodic boundary condition ([Wisdom &: Tremaind Il988l ). There are an active 
region and its copies. Inner and outer regions have different angular velocities. These regions 
slide upward and downward with shear velocity 3r2o-^x/2, where Lx is the length of the cell in 
the x-direction. We calculate gravitational forces from particles in these regions. We use the 
special-purpose computer GRAPE-6 for the calculation of self-gravity (iMakino et al.l 120031 ). 



Particles used in the simulation are not realistic dust particles but superparticles that represent 
a group of small particles, thus the parameter e used in this Paper is not exactly the physical 
restitution coefficient but corresponds to the coefficient of the dissipation due to collisions among 
particles. The collision models used in paper I cannot treat more dissipative models where the 
tangential coefficients of restitution et < 1 because we adopt et = 1. The net restitution coefficient 
is e = (en(fn) -|- f-t {vt)) / {{Vn) + (^t))) where is the normal coefficient of restitution, and ft 
are normal and tangential velocities of a collision (ICanup &: Espositdll995l ). The net restitution 
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coefficient e is always larger than 0.7 if we assume et = 1 and (w^) = (ft)- If we allow et 7^ 1 in 
order to examine the model of e < 0.7, we must consider rotations of individual particles. For the 
sake of simplicity, we ignore rotations of particles in this paper and assume et = 1. 

The way to treat a collision is to use the same method applied to the hard sphere model, but 
we consider the binding condition of colliding particles. If a pair overlaps, ftj < fp,i + fpj where fp^j 
is a radius of particle i, and is approaching, fij,- • Vj,- < where Vj,- is the relative velocity and fij,- 



is a unit vector along f = {xj , yj ,Zj) — {xi ,yi,Zi 
its velocities by using the equations of collision: 



we assume that this pair collides, and changes 



rrii 



rhn + m. 



■(1 + e){i\ij ■ Vij)n: 



(2) 



Vi + - 



rrii 



(3) 



where e is a restitution coefficient. Using the u 
condition is satisfied by the Jacobi integral (e.g.. 



J, 



2 



+ Vij + z, 



~ 2. 



jdated velocities , we c heck whether the binding 
+ ^ < 0, (4) 



Nakazawa fc Idalligsa ): 
3 



3 _9 1 _9 
— T- • -I 

2 V 2 



where {xijjy^j, Zij) is the relative velocity. If the pair is bound, two particles are merged into one 
particle conserving their momentum. The shape of the new particle is a sphere, the radius of which 
is determined by its mass. 

We assume that all particles have the same mass rrip initially. The parameters of the present 
simulation are the distance from the central star oq, the length of region Lx and Ly, the number 
of particles N, and the mass of particles mp(rp,/9p), where pp is the densi ty of a particle. Th e 
dynamical behavior is characterized by only two non-dimensional parameters ([Paisaka &: Idalll999l ): 
the initial optical depth r and the ratio ( = aoh/2rp. In the models of ^ > Ij the centers of the 
particles are i n their Hill sphere when two particles come into contact. The initial optical depth is 
given by (e.g.. iGoldreich &: Tremaindll982l ) 



3S 



0.19 



4/9prp \10gcm \2gcm •^y V20cm 

where S is the surface density of particles. The ratio C is given by 



(5) 



C = 105.528 



Mo 



^1/3 



Pp 



2gcm 



-3 



1/3 



lAU/ ■ 



(6) 



Using the above two parameters, the other paramet ers can be es timated. The normalized most 
unstable wavelength of gravitational instability is (e.g.. ISekiyalll983l ) 



Amost = 127rTC^. 



(7) 
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The normalized size of computational domain L, the number of particles N, the normalized radius 
of particle fp, and Toomre's Q value are given by 



= UttA^tC"^ 



N 



3/-6 



576TTA^AyT-^C 

^ _ 1 

"^P" 2C' 
Q 



6tC 



(8) 

(9) 
(10) 

(11) 
(12) 



where A^ is the ratio Lx/Xmost, Ay is the ratio Ly /X^ost, and ax is the radial velocity dispersion 
scaled by aoh^lo, respectively. In order to determine the size of the computational domain, we need 
to set Ax and Ay. The initial velocity dispersion is also a parameter. We determine it from the 
initial Toomre's Qinit value using Equation (fT2]) . As shown above, if we set r, Ax, Ay, and Qinit, 
we can determine the initial state. The simulation parameters are summarized in Table [TJ We call 
model 1 as the standard model. The number of particles is 1000-10000. 



If C > 1.5, two particles can be bound (|Qhtsukilll993l : ISaldll995l ). so we need to adopt ( > 1.5 
in order to investigate the formation process of planetesimals. In this paper, we use C = 2, which is 
much smaller than the realistic value ( ~ 100. Thus the physical size of a particle in our simulation 
is very large. It indicates that the effect of the self-gravity is relatively weaker than that of inelastic 
collisions. To investigate gravitational instability, the particle density should be higher than the 
Roche density. The particle density is pp = A^^p^i, where pR = (9/47r)M*/aQ is the Roche density. 
Thus when = 2, pp ^ /jr. So we can study gravitational instability with this parameter. But 
this approximation may change the result when comparing between C = 2 and ^ = 100 although 
gravitational instability occurs in both models. In the future work, we will perform the simulation 
with a larger C value by using next-generation supercomputers. 

In most of the models, we set Qinit = 3, which corresponds to ax — 7.2 for the standard model. 
The characteristic velocity of the actual turbulence fturb is about tjvk, where vk is the Kepler 
velocity, rj = — 12(c^/w^) (Slog P/5 log r), c,, is the sound velocity of gas, P is the p ressure of gas. 



tne_pr( 
197^ ). 



and r is the distance from the central star (jAdachi et al.lll976l : IWeidenschillinall977l ). At r = lAU, 
in the Hill unit, the characteristic velocity of the turbulence is Uturb — 1-3 x 10'' that is much 
larger than ax- If the turbulence weakens, Q decreases gradually from Q ^ 1, and gravitational 
instability occurs finally. Therefore we choose the initial condition that the disk is gravitationally 
stable (Qinit = 3). 

We set the initial velocities and positions of particles using the above parameters. A Keplerian 
orbit has six parameters: the position of the guiding center (xg,?/g), eccentr icity e, inclination i, and 
the two phases of epicycle and vertical motions (e.g.. lNakazawa &: Idalll988l ). The eccentricity e and 
inclination i of particles are assumed to obey the Rayleigh distribution with the ratio {e^) / (z^) = 2 
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Table 1. Initial conditions and Results 



Model 


e 


r 


c 




Ay 


Qinit 




(Mpi/Miinear) 


iVpi 


Ul 


U.Ul 


u.l 


^.U 


O.U 


O.U 


O.U 


7.Z 


21. z 


2 


no 


U.Ul 


0.11 




O.U 


O.U 


O.U 


T A 

7.9 


4:2A 


1 


no 


U.Ul 


n 1 occ 
O.lzo 


z.U 


O.U 


O.U 


O.U 


y.u 


21.0 


o 
2 


(J4 


U.l 


U.l 


z.U 


b.U 


b.U 


3.U 


7.2 


42.7 


1 


Do 


U.3 


U.l 


z.U 


b.U 


b.U 


3.U 


7.2 


/loo 


1 


Ud 


U.5 


U.l 


2.U 


O.U 


b.U 


3.U 


7.Z 


o8.4 


1 


U7 


U.7 


U.l 


2.U 


O.U 


O.U 


O.U 


7.Z 


A A 

u.u 


U 


no 
Do 


U.9 


U.l 


Z.U 


O.U 


b.U 


O.U 


7.Z 


A A 

U.U 


U 


uy 


U.Ul 


U.Ub 


2.0 


O.U 


b.U 


O.U 


b.8 


AO A 

4o.4 


1 


lU 


U.Ul 


U.Ub 


Z.I 


b.U 


b.U 


O.U 


8.2 


AO A 

4J.4 


1 


ii 


U.Ui 


U.Ub 


3.U 


b.U 


b.U 


O A 

O.U 


9.7 


OA O 

o9.8 


1 


12 


0.01 


U.l 


2.U 


4.U 


4.U 


3.U 


7.2 


19.5 


1 


13 


0.01 


0.1 


2.0 


8.0 


8.0 


3.0 


7.2 


38.4 


2 


14 


0.01 


0.1 


2.0 


3.0 


3.0 


3.0 


7.2 


11.1 


1 


15 


0.01 


0.1 


2.0 


6.0 


3.0 


3.0 


7.2 


21.7 


1 


16 


0.01 


0.1 


2.0 


9.0 


3.0 


3.0 


7.2 


16.4 


2 


17 


0.01 


0.1 


2.0 


12.0 


3.0 


3.0 


7.2 


11.0 


4 


18 


0.01 


0.1 


2.0 


24.0 


3.0 


3.0 


7.2 


10.1 


9 


19 


0.01 


0.1 


2.0 


3.0 


6.0 


3.0 


7.2 


21.8 


1 


20 


0.01 


0.1 


2.0 


3.0 


9.0 


3.0 


7.2 


32.2 


1 


21 


0.01 


0.1 


2.0 


3.0 


12.0 


3.0 


7.2 


43.1 


1 


22 


0.01 


0.1 


2.0 


3.0 


24.0 


3.0 


7.2 


85.5 


1 


23 


0.01 


0.1 


2.0 


6.0 


6.0 


4.0 


9.6 


43.3 


1 


24 


0.01 


0.1 


2.0 


6.0 


6.0 


5.0 


12.0 


21.6 


2 



Note. — (-^pi)) -^linear and Npi are the mean planetesimal mass, the linear 
mass and the number of planetesimals at final state, respectively. 
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(|lda k Makinol 11992 1 
avoiding overlapping. 



We set the position of the guiding center and the two phases uniformly, 



The equation of motion is integrated w ith a second-order leapfrog scheme. We adopt the vari- 



able time step used by lDaisaka &: Idal (119991 ). The time step formula is given by At = miuj |aj|/|aj 



where aj and a^ are the acceleration of particle i and its time derivative, and rj is a non-dimensional 
time step coefficient, respectively. 



3. Results 

3.1. Time Evolution 

Figure [I shows snapshots in model 1 at t/Tx = 0.0,0.4,0.8,1.2,1.6 and 2.0 where Tk is the 
orbital period 2tt/Qq. At t/Tx = 0.4, we cannot see any spatial structures and large bodies. The 
formation process of planetesimals is the same as those of the hard and soft sphere models (Paper 
I). The non-axisymmetric wake-like structure forms at t/Tx = 0.8, 1.2, which appears if we conside r 
self-gravity. The non-axisymmetric wake-like structure is also seen in planetary rings (|Salol[l995l ). 
The density in the wake-like structure is higher than the mean surface density S. For example, the 
surface density in the dense region is 1.52S at t/Tx = 1.6. Planetesimal seeds form in the dense 
parts of these wakes by fragmentation (t/Tx = 1.6,2.0). Here we consider the formation of the 
wake-like structure and planetesimal seeds as the gravitational instability stage. 

Figure [2] shows snapshots in model 1 at t/Tx = 2.0, 3.0, 4.0, 5.0, 6.0 and 7.0. Once planetesimal 
seeds form, the non-axisymmetric wake-like structures disappears (t/Tx = 3.0). Planetesimal seeds 
grow rapidly by mutual coalescence (t/Tx = 3.0,4.0,5.0). In this stage, the disk is gravitationally 
stable. The gravitation seems to merely enhance the coalescent growth. When almost all plan- 
etesimal seeds are absorbed by a few planetesimals, the growth slows down (t/Tx = 6.0, 7.0). This 
stage is the collisional growth stage. 

Figure [3] shows the time evolution of the number of large particles, the mass of the largest and 
second largest particles, the ratio of the mass of large particles to the total mass and the velocity 
dispersion of field particles for model 1. The time evolution of the number of large particles is 
similar to those of the hard and soft sphere models (Paper I). The number of large particles has a 
maximum at t/Tx — 1.7, and its value is about 1.2, which is slightly larger than that in the hard 
and soft sphere models. However, the number of large particles decreases rapidly after t/Tx — 1.7. 
This decay is caused by the fast coalescence of large particles. 

The mass of the largest particle monotonically increases up to about 44Miinear; which is twice 
larger than that of the hard and soft sphere models, where Miinear = 7rS(Ainost/2)^ is the plan- 
etesimal mass estimated by the linear theory. In the accretion model, the growth of planetesimals 
is enhanced because sticking of dust grains is efficient. The mass of the second largest particle 
changes discontinuously, when the second largest particle is absorbed by the largest one. 
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The ratio of the mass of the large particles to the total mass is almost the same as that of the 
hard and soft sphere models. The ratio monotonically increases up to about 1, and most of the 
mass in the computational domain is finally absorbed by the large planetesimals. 

The velocity dispersion decreases until t/Tx ~ 1.3 because of the dissipation due to inelastic 
collisions. When the velocity dispersion becomes sufficiently small {Q ~ 2), the gravitational 
instability occurs and a lot of planetesimal seeds form. As the field particles are scattered by the 
planetesimal seeds, the velocity dispersion increases. 



The strength of self-gravity is measured by two dimensionless quantities: (e.g.. lYoudinll2005al ) 



Q^^. (13) 

where p is the mass density of the dust layer, and is the thickness of the dust layer. The value 
Qr is the ratio of the Roche density to the dust density. If these values are sufficiently small, 
the disk is gravitationally unstable. Figure [4] shows time evolution of Q and Qr. The value Qr 
is always larger than Toomre's Q value. The ratio Q/Qr is approximately equal to 0.5. If the 
hydrostatic balance {ax — ^ohd) can be assumed, Q/Qr = 4/9 ~ 0.44. 

If the gravitational instability occurs, there are no remarkable qualitative differences in the 
planetesimal formation process among different collision models. The growth of particles in the 
accretion model is more efficient than in the hard and soft sphere models, thus the maximum mass 
of planetesimals is slightly larger than that in the hard and soft sphere models. 



3.2. Parameter Dependence 

3.2.1. Physical Parameters 



We summarize the parameter dependences of the simulation results on the restitution coeffi- 
cient e, the optical depth r, the ratio of Hill radius to the particle diameter C,, the initial Toomre's 
Q value. We focus on the mass of the largest particle, and Toomre's Q value. 

Figure [5] shows the dependence on the restitution coefficient e. We vary the restitution coef- 
ficient e from 0.01 to 0.9 (models 1,4,5,6,7,8). No gravitational instability occurs, in other words, 
Q never becomes less than about 2, for e = 0.9, and no large body form. There are two reasons 
why planetesimals do not f orm. The velocity dispersion i ncreases rnonotonically if e is la r ger than 
the c ritical value ^ 0.7 (jColdreich k Tremainj|l982l : ISalolll995l : IPaisaka fc Idalll999l : lohtsuki 
19991 ). The velocity dispersion decreases due to inelastic collisions and increases due to gravita- 



tional scattering. If the restitution coefficient is large, the dissipation is relatively weak, thus the 
velocity dispersion does not decrease. Because Q value does not reach the critical value, gravita- 
tional instability does not occur. Therefore, planetesimal seeds do not form, and the subsequent 
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Fig. 1. — Spatial distribution of particles in the standard model (model 1) at t/Tx = 0.0 {top left 

panel), t/Tn = 0.4 {top middle panel), t/T^ = 0.8 {top right panel), t/T^ = 1.2 {bottom left panel), 
t/Tiii = 1.6 {bottom middle panel), and t/Tx = 2.0 {bottom right panel). Circles represent particles 
and their size is proportional to the physical size of particles. 
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Fig. 2. — Same as Figure[I]at t/Tx = 2.0 {top left panel), I/Tk = 3.0 {top middle panel), t/T^ = 4.0 
{top right panel), t/T^ = 5.0 {bottom left panel), t/T]^ = 6.0 {bottom middle panel), and t/Tx = 7.0 
{bottom right panel). 
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Fig. 3. — Time evolution of the various quantities in the standard model (model 1). The quantities 
are the number of particles with m > O-lMunear per the area A^^gt, the mass of the largest and 
second largest bodies Mmax normalized by the linear mass Muenar, the ratio of the total mass of 
particles with m > OTMnnear to total mass Mtotai and the velocity dispersion of field particles from 
top to bottom, respectively. In the bottom panel, we show the velocities for Q = 2 (dotted line) 
and Q = 1 {dashed line) 
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collisional growth does not happen. We cannot apply the formation process of planetesimal stated 
in our paper to the model with e = 0.9. With a larger restitution coefficient, the reduction of the 
relative velocity on collisions is smaller. In addition, the collision velocity increases as the velocity 
dispersion increases. The coagulation is not efficient with the large restitution coefficient. Thus 
planetesimals do not form by coagulation within the simulation time. 

We study the effect of the optical depth r by comparing results with r = 0.1,0.11, and 0.125 
(models 1,2,3). The dependence on the optical depth r is shown in Figure[6l The results are similar 
to those of the hard and soft sphere models (Paper I). The largest masses and Toomre's Q value 
are on the same order. 

We set C, = 2.5, 2.75, and 3.0 and fix the other parameters in the standard model except for 
the optical depth r = 0.06 (models 9,10,11). The dependence on Q is shown in Figure El . The 
result is also similar to that of the hard and soft sphere models (Paper I). The largest masses are 
A^max /-^linear = 43.3,43.2, and 39.8 respectively. The largest masses are almost the same, which is 
about 40Miinear- Figure [7] shows the similar time evolution of Toomre's Q value because the masses 
of planetesimals are almost the same. 

Figure [8] shows the time evolution of Toomre's Q values for Qmit = 3.0, 4.0, and 5.0 (models 
1, 23, 24). If the growth is due to gravitational instability, the results do not depend on the initial 
Toomre's Q value Qinit- Time evolutions are similar for these models. Toomre's Q value decreases 
through inelastic collisions between particles. As Q decreases, self-gravity becomes stronger. When 
Q becomes less than about 2, the disk becomes gravitationally unstable and the wake-like structure 
develops. Then Q increases due to scattering of field particles by planetesimal seeds. 

In Paper I, we showed that there is no remarkable difference between the hard and soft sphere 
models on the formation process of planetesimals. We confirm that the results of the accretion 
model are also essentially the same as these of the hard and soft sphere models. The dissipation 
rate due to collision is a crucial parameter but the planetesimal formation process is independent 
of the collision model. Before the gravitational instability, collisions merely decrease the velocity 
dispersion. In the late stage, coalescence among planetesimal seeds occurs, but this process is 
independent of collision models. 

3.2.2. Simulation Domain Size 

In the hard and soft sphere models, the mass of the largest planetesimal depends on the size of 
the computational domain, and we could not find the saturation of growth (Paper I). The growth 
of planetesimals continues until most of the mass in the computational domain is absorbed by 
the largest planetesimal. We expect to find the saturation of growth because the accretion model 
enables us to perform a wider simulation with a large number of particles. 

First we assume that the ratio of the length of the computational domain in the x-direction 
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Fig. 4. — Time evolution of Toomre's Q value Q (solid line) and the ratio of the Roche density to 
the dust density Qr (dashed line) for model 1. 
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Fig. 5. — Time evolution of various quantities for e = 0.01 {solid curve), e = 0.1 {dashed curve), 
€ = 0.3 {short dashed curve) , e = 0.5 {dotted curve), e = 0.7 {dash-dotted curve), e = 0.9 {dot-short- 
dashed curve) (models 1,4,5,6,7,8). The quantities are the mass of the largest body, and Toomre's 
Q value from top to bottom, respectively. 
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Fig. 6. — Same as Figure [5] but for r = 0.1 {solid curve), t = 0.11 (dashed curve), t = 0.125 
{dotted curve) (models 1,2,3). 




Fig. 7. — Same as Figure [5] but for ( = 2.5 {solid curve), C = 2.75 {dashed curve), ( = 3.0 {dotted 
curve) (models 9,10,11). 
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Fig. 8. — Time evolution of Toomrc's Q value. Initial Q values are 3.0 {solid curve), 4.0 {dashed 
curve), 5.0 {shot dashed curve) (models 1, 23, 24). 
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to the y-direction is unity, in other words, the shape of the computational domain is a square. We 
vary the size of the computational domain ^ = 4 to 8. Results are shown in Figure [9] (models 12, 
1, 13). The final masses of the largest planetesimals and Toomre's Q values are Mmax /-^linear = 
19.5,44.0,66.5, Q = 13.6,18.2,21.9, respectively. As the size of the simulation domain increases, 
the final planetesimal mass and the final Toomre's Q value increase. 

Next, we vary the ratio of the length of the computational domain in the x- and y-directions. 
Now the shape of the computational domain is rectangular. We vary the length in the y-direction. 
Ay, from 3 to 24 fixing the length in the x-direction. Ax, (the left panel in Figure [101) (mod- 
els 14, 15, 16, 17, 18), and we vary A^ from 3 to 24 fixing Ay (the right panel in Figure [TO]) 
(models 14, 19, 20, 21, 22). The masses of the largest planetesimal and Toomre's Q values are 
Mmax/Miinear = H-l, 21.7, 18.8, 15.9, 16.8, Q = 10.4, 13.6, 10.2, 12.2, 11.3 (models 14, 15, 16, 17, 18) 
and Mmax/Miinear = 11.1,21.8,32.2,43.1,85.5, Q = 10.4,17.1,15.9,21.3,24.1 (models 14, 19, 20, 
21, 22), respectively. In the model with fixed Ax , the mass of the largest planetesimal increases 
with Ay, while in the model with fixed Ay, the largest mass is roughly constant. These results 
indicate that the mass of the largest planetesimal is determined by Ay. If Ax is sufficiently large, 
multiple planetesimals form. Toomre's Q value increases with Ay because of the scattering of field 
particles by the large bodies. 

Here we estimate the final mass of a planetesimal, Mpi, in the computational don iain in a similar 
way t o estimate the isolation mass of protoplanets in the oligarchic growth picture (iKokubo &: Ida 



19981 ). We assume that particles whose orbital radii are within of the planetesimal are finally 
absorbed by the planetesimal though they are not in the Hill sphere of the planetesimal initially, 
where 7 is a non-dimensional factor and rn is the Hill radius of the planetesimal. Therefore, the 
planetesimal absorbs particles in the rectangle of 27rH in width and Ly in length. The hill radius 
of the planetesimal is rn = ao(2Afpi/3M,)i/3, and we define 

Mh = 2-trnLyT.. (15) 




Fig. 9. — Same as Figure [5] but for Ax = Ay = A {solid curve), Ax = Ay = Q [dashed curve), and 
Ax = Ay = S {dash-dotted curve) (models 12,1,13) . 
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The planetesimal with the mass Mpi can absorb the mass Mh , and Mh increases as the planetesimal 
grows. When Mpi ~ Mjj, the growth becomes slow. We estimate the mass of the planetesimal by 
this condition. In this case the number of planetesimals is about Lx/ {2'yr}i). We obtained the mass 
of the planetesimal 

Mpi = min (bnV^er\^Ayfl\M,,,^^ , (16) 



where Mtotai are the total mass in the computational domain. The number of the planetesimals is 
A^pi is 

iVpi = max(^vry|-^A^,l). (17) 

The mass of the planetesimals is equal to the total mass in the computational domain and the 
number of the planetesimal is unity if the following condition is satisfied: 



A. < ^<%^/^. (18) 



Svr 

The comparisons between the analytical estimation and the simulations are shown in Figures 
[TT] and [T2j If we assume 7 = 2.5, the analytical estimation agrees with the simulation results. 
The parameter 7 = 2.5 indicates that the mean orbital separation of planetesimals is 5rH. This 
is narro wer than lOrH, which is the typical orbital separation of protoplanets caused by oligarchic 



growth (iKokubo &: Idalll998l . 12000) • This narrow orbital separation is possible for a dissipative 
system. We will discuss this issue later. For instance, Figure [13] shows the final state in model 
18. In this Figure, there are 9 planetesimals in the computational domain whose length in the 
x-direction is 400. Thus the mean orbital separation is about 44. The Hill radius of the mean 
planetesimal ((Mpi/Miinear) — 10.1) is about 9.5. Therefore the separation is estimated to about 
9.5 X 5 = 48, which is consistent with Figure [T3l 

Our analytical estimation gives good agreement with the numerical results in each model. 
Figure [TT] shows the mean mass and the number of planetesimals as a function of the size of the 
computational domain in the model where the domain is square. The mean mass of planetesimals 
increases with the size of the computational domain. From the Equation (jl8p . when = Ay < 4.2, 
the number of planetesimals is equal to unity and the mass of the planetesimal is approximately 
equal to the total mass in the computational domain. The slope of the line of the mean mass of 
planetesimals changes at Ax = Ay = 4.2 and the number of the planetesimals is larger than unity 
if Ax = Ay> 4.2. 

Figure [12] shows the result for the model of the rectangular computational domain. In the 
models of Ay = 3, the mean mass of planetesimals is constant and the number of planetesimals is 
larger than unity if A^ > 3.6. In the models of Ax = 3, the slopes of the lines change at Ay = 2.1. 
The growth of planetesimals cannot be saturated although Ay is large. In these models, the Hill 
radius of the planetesimal is longer than the size of the computational domain in the x-direction. 
The only one planetesimal forms, and it absorbs most mass in the computational domain. It 
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should be remembered that the local calculation of planetesimal formation has the domain-size 
dependence. 



If we assume that the surface density of dust is 

S 

where fd is a scaling factor, the linear mass Mnnear is given by 



10 X /d ( ) gem , 



Miinear = 3.5 X 10 
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and from Equation (jl6p the estimated planetesimal mass Mp\ is given by 



Mpi = 9.0 X lO^Vd 
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where we assume > 2.1AI'^ (7/2.5)^/^ Th e factor /h = 1 roughly corresponds to the minimum 



mass solar nebula model inside the snow line (|Hayashilll98ll ). From Equation (I2ip . the estimated 
planetesimal mass becomes larger than the linear mass if the feeding zone is large, (Ay'y)^^'^ > 1.6. 

The estimated and linear masses exhibit the same dependencies on /d, ao, and M^:, and the 
estimated planetesimal mass is on the same order of the linear mass when Ay'j ~ 1. These results 
are explained as follows. Here we rewrite the Hill radius by the most unstable wavelength: 
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most 



( Afpl \ 



1/3 
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where the most unstable wavelength for Q = 1 is 

27r2GE 
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(22) 
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The most unstable wavelength of the gravitational instability and the Hill radius of the linear mass 
are on the same order. The condition Q = 1 corresponds to the balance between the self-gravity 
and the effect of the rotation. On the other hand, the Hill radius is determined by the balance 
between the self-gravity and the tidal force. Thus, rn and Amost are on the same order, and the 
planetesimal mass is on the order of the linear mass only if Ayj ~ 1. In general, because Ay ^ 1, 
the planetesimal mass is larger than the linear mass. 

Collisions among planetesimals are more frequent than the realistic case in the present simu- 
lations since the bulk density of super-particle is much smaller than the realistic value. Therefore 
the orbital separation of large planetesimals is di fferent from that of the sta ndard oligarchic growth 
derived by A^-body simulations of planetesimals (jKokubo &: Idalll998l . l2000l ). If the number density 
is low, t he orbital separation is determined by the balance between the orbital repulsion and their 
growth (jKokubo Idalll998l ). But if the number density is high, and the collisional dissipation 
is effective, the orbital repulsion is relatively w eak. Thus, in the coll isional oligarchic growth, the 
orbital separation can be narrower than IORh (jGoldreich et al.ll2004l ). 
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Fig. 10. — Same as Figure [5] but for Ax = 3 {solid curve), = G {dashed curve), A^ = 9 {dash- 
dotted curve). Ax = 12 {short dashed curve), and Ax = 24 {dotted curve) (models 14,15,16,17,18) 
in the left panel, and Ay = 3 {solid curve). Ay = 6 {dashed curve). Ay = 9 {dash-dotted curve). 
Ay = 12 {short dashed curve), and Ay = 24 {dotted curve) (models 14,19,20,21,22) in the right 
panel. 
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Fig. 11. — Averaged mass {Mpi) {left panel) and the number A^pi {right panel) of the planetesimals 
as a function of the size of the computational domain A = Ax = Ay = 4:,6,8 (models 12,1,13). The 

2.5. The 



solid line denotes the analytical estimation given by Equations (|T6ll and (|T7|) . We set 7 
error bar shows the range of planetesimal masses in a given run. 
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Fig. 12.— Same as Figure E] but for A = = 3, 6, 9, 12, 24 {cross points) (models 14,15,16,17,18) 
and A = Ay = 3,6, 9, 12, 24 {plus points) (models 14,19,20,21,22). The dotted and solid lines denote 
the analytical estimations given by Equations (|16p and (jl7p . 
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Fig. 13. — Spatial distribution of particles in model 18 at t = AOTk- The mean planetesimal mass 
is about (Mpi/Miincar) — 10 a-iid its Hill radius is about 9.5. 
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4. Summary 

We performed local A^-body simulations of planetesimal formation through gravitational insta- 
bility and collisional growth without a gas component using a shearing box method. The accretion 
model was adopted as the collision model. In the accretion model, when particles collide, if the 
binding condition is satisfied, two particles merge into one. The number of particles decreases as 
the calculation proceeds, thus this model enables us to perform the long-term and large-scale cal- 
culations. We compared the results of the accretion model with those of the hard and soft sphere 
models used in Paper I. 

The formation process of planetesimals is the same as that of the hard and soft sphere mod- 
els. It is divided into three stages: the formation of non-axisymmetric structures, the creation 
of planetesimal seeds, and the collisional growth of the planetesimals. The velocity dispersion of 
particles decreases gradually due to the inelastic collisions. The gravitational instability occurs 
when the velocity dispersion of particles becomes a critical value determined by Q ~ 2. The non- 
axisymmetric structures form and fragment into planetesimal seeds, the masses of which are on the 
order of the linear mass. In this paper, we consider the formation of non-axisymmetric structures 
and planetesimal seeds as gravitational instability. Planetesimal seeds grow rapidly by coagulation 
and the large planetesimals form. Most of the mass in the computational domain is absorbed by a 
few large planetesimals. 

We studied the effect of physical parameters on planetesimal formation: the restitution co- 
efficient e, the initial optical depth r, the dependence of planetesimal formation on the ratio 
= rH/2rp, and the initial Toomre's Q value. We found no qualitative differences among the 
collision models. In the accretion model, the merged particles cannot split, thus the growth of 
particles is more efficient than in the hard and soft sphere models. If the restitution coefficient e is 
smaller than the critical value, the time evolution is similar to the standard model. If r is large, the 
collision frequency is high, and the dissipation of the kinetic energy is efficient. Thus, the gravita- 
tional instability occurs more quickly for larger r. In this narrow parameter range ( = 2.5, 2.75, 3.0, 
there is no remarkable dependence on 

By the long-term and large-scale calculations using the accretion model, we found that the 
mean mass of the planetesimals depends on the size of the computational domain. We found that 

3 /2 

the mean mass of planetesimals is proportional to Ly . However, this mass is independent of 
if Lx is sufficiently large. The mean mass of planetesimals is estimated by Equation (jl6p . Large 
planetesimals sweep small planetesimals in the rotational direction. If the orbital separation of 
planetesimals is sufficiently large, the planetesimals cannot collide. The typical orbital separation 
is about 5rH in the present simulation The dependence of the planetesimal mass on the size of the 
computational domain indicates that we cannot simply discuss the realistic mass of the planetesimal 
using the local calculation. 

The effect of gas was neglected in our simulation in order to study the physical process of 
gravitational instability of dust particles and subsequent collisional evolution as a first step. Once 
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gravitational instability happens and large aggregates form, we may be able to neglect gas drag 
because they are so large that they decouple from gas. So our gas-free model may be applicable 
to the stages of gravitational instability and subsequent collisional growth. The simple gas-free 
model provides thorough understanding of the gravitational instability and collisional growth of 
planetesimals. However, it is necessary to include gas drag so as to investigate the realistic formation 
process of planetesimals since the stopping time of small dust particles is shorter than the Kepler 
time. We will investigate the gravitational instability of a particle disk embedded in a laminar 
gaseous disk in the next paper. 
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